1 Introduction

Linkage maps are essential tools to several genetic endeavours such as quantitative trait loci (QTL) analysis, evolutionary studies, assessment of colinearity in chromosome-wise scale between genomes, and study of meiotic processes. The principle behind linkage map construction is detecting recombinant events between genomic positions and summarizing them into pairwise recombination fraction estimates. In diploids, the assessment of such a phenomenon is relatively straightforward. After the homolog duplication, sister chromatids pair up exchanging segments. The presence of informative markers, e. g. single nucleotide polymorphisms (SNPs) enable computation of the recombination fraction between pairs of genomic positions by comparing the chromosome constitution of parents and offspring. If the order of markers is unknown, it can be obtained using the pairwise recombination fractions in conjunction with optimization algorithms.

In polyploids (species with more than two sets of homolog chromosomes, or homologs), the construction of such maps is quite challenging. While in diploids a biallelic marker can be fully informative, in polyploids they allow accounting for proportions of biallelic dosages. To recover the multiallelic information present in polyploid species, we need to account for recombination frequencies, estimate phase configurations and reconstruct the haplotypes of both parents and individuals in the population. As the ploidy level increases, the joint computation of genomic regions becomes intensive, and other approaches such as dimension reduction need to be applied.

MAPpoly is an R package fully capable of building genetic linkage maps for biparental populations in polyploid species with ploidy level up to 8x. This package is part of the Genomic Tools for Sweetpotato Improvement (GT4SP), funded by Bill and Melinda Gates Foundation. All of the procedures presented in this document are detailed in Mollinari and Garcia (2019) and Mollinari et al. (2020). The results obtained with MAPpoly can be readily used for QTL mapping with the QTLpoly package, which implements the procedures proposed by da Silva Pereira et al. (2020).

Some advantages of MAPpoly are:

• Can handle multiple dataset types
• Can handle thousands of markers
• Does not depend on single-dose markers (SDM) to build the map
• Incorporates genomic information
• Explores multipoint information (through Hidden Markov Chain Model)
• Can handle high ploidies: up to 8x when using HMM, and up to 12x when using the two-point approach
• Can reconstruct haplotypes for parents and all individuals in the population
• Recovers the multiallelic nature of polyploid genomes
• Detects occurrence, location and frequency of multivalent pairing during meiosis
• Robust enough to build genetic linkage maps with multiallelic markers (when available)

install.packages("mappoly")

install.packages("devtools")

devtools::install_github("mmollina/mappoly", dependencies=TRUE)

2 Loading datasets

In its current version, MAPpoly can handle the following types of datasets:

1. CSV files
2. MAPpoly files
   • Dosage based
   • Probability based
3. fitPoly files
4. VCF files

MAPpoly also is capable of importing objects generated by the following R packages

1. updog
2. polyRAD
3. polymapR
   • Datasets
   • Maps

Both CSV and MAPpoly datasets are sensible to formatting errors, such as additional spaces, commas and wrong encoding (non-UTF-8). If you have any trouble, please double check your files before submitting an issue. You can find detailed steps of all supported files in the following sections. Also, in large datasets, it is expected that a considerable proportion of markers will contain redundant information. Thus, all reading functions will set these redundant markers aside to be included in the final map.

2.1 Reading CSV files

The preparation of a CSV file for MAPpoly is relatively straightforward. It can be done in Microsoft Excel or any other spreadsheet software of your preference. In this file, each line comprehends a marker and each column comprehends information about the marker. In its current version, MAPpoly can handle .csv files with allelic dosage data.

The first line of the CSV file should contain headers for all columns. The first five columns should include the following information: marker name, dosage for both parents (one column for each), a sequence number (e.g. a chromosome number, if available) and a sequence position (e.g. the marker position within the chromosome, if available). In addition to these five headers, you should include the name of all individuals in the population. From the second line onwards, all columns should contain its values, including allelic dosages for all individuals. Missing or absent values should be represented by NA.

NOTE: If genomic information is not available, the ‘sequence’ and ‘sequence position’ columns should be filled with NA’s.

Example:

Figure 1: Example of CSV data set

Important note: avoid spaces in .csv files. As mentioned above, please double check your datasets for extra spaces, commas, dots and encoding. Your CSV file should be encoded using UTF-8.

You can read CSV files with the read_geno_csv function:

ft="https://raw.githubusercontent.com/mmollina/SCRI/main/data/B2721_dose.csv"
tempfl <- tempfile()
download.file(ft, destfile = tempfl)
dat.dose.csv <- read_geno_csv(file.in  = tempfl, ploidy = 4)
#> Reading the following data:
#>     Ploidy level: 4
#>     No. individuals:  156
#>     No. markers:  7093
#>     No. informative markers:  7093 (100%)
#>     This dataset contains sequence information.
#>     ...
#>     Done with reading.
#>     Filtering non-conforming markers.
#>     ...
#>     Done with filtering.
unlink(tempfl)

dat.dose.csv

In addition to the CSV file path, you must indicate the ploidy level using the ploidy argument. This function automatically excludes uninformative markers. It also performs chi-square tests for all markers, considering the expected segregation under Mendelian inheritance, random chromosome pairing and no double reduction. You can optionally use the filter.non.conforming logical argument (default = TRUE), which excludes non-expected genotypes under these assumptions. However, keep in mind that the linkage analysis functions in MAPpoly do not support double-reduction in its current version, and this can cause software failures.

ploidy 4
nind 3
nmrk 5
mrknames M1 M2 M3 M4 M5
indnames Ind1 Ind2 Ind3
dosageP 0 2 0 0 3
dosageQ 1 2 1 1 3
seq 1 1 2 2 3
seqpos 100 200 50 150 80
nphen 0
pheno-----------------------
geno------------------------

ploidy 4
nind 3
nmrk 5
mrknames M1 M2 M3 M4 M5
indnames Ind1 Ind2 Ind3
dosageP 0 2 0 0 3
dosageQ 1 2 1 1 3
seq 1 1 2 2 3
seqpos 100 200 50 150 80
nphen 0
pheno-----------------------
geno------------------------
1 0 0
3 0 2
1 0 0
1 0 0
3 4 4

fl = "https://raw.githubusercontent.com/mmollina/SCRI/main/data/B2721_mappoly_dose"
tempfl <- tempfile()
download.file(fl, destfile = tempfl)
dat.dose.mpl <- read_geno(file.in  = tempfl, elim.redundant = TRUE)
#> Reading the following data:
#>     Ploidy level: 4
#>     No. individuals:  156
#>     No. markers:  6502
#>     No. informative markers:  6502 (100%)
#>     This dataset contains sequence information.
#>     ...
#> 
#>     Done with reading.
#>     Filtering non-conforming markers.
#>     ...
#>     Performing chi-square test.
#>     ...
#>     Done.
unlink(tempfl)

dat.dose.mpl

ploidy 4
nind 3
nmrk 5
mrknames M1 M2 M3 M4 M5
indnames Ind1 Ind2 Ind3
dosageP 0 2 0 0 3
dosageQ 1 2 1 1 3
seq 1 1 2 2 3
seqpos 100 200 50 150 80
nphen 0
pheno-----------------------
geno------------------------
M1 Ind1 0.5 0.5 0.0 0.0 0.0
M2 Ind1 0.0 1.0 0.0 0.0 0.0
M3 Ind1 0.3 0.7 0.0 0.0 0.0
M4 Ind1 0.5 0.5 0.0 0.0 0.0
M5 Ind1 0.0 0.0 0.0 0.9 0.1
M1 Ind2 1.0 0.0 0.0 0.0 0.0
M2 Ind2 0.2 0.5 0.3 0.0 0.0
M3 Ind2 0.9 0.1 0.0 0.0 0.0
M4 Ind2 0.9 0.1 0.0 0.0 0.0
M5 Ind2 0.0 0.0 0.0 0.2 0.8
M1 Ind3 0.2 0.8 0.0 0.0 0.0
M2 Ind3 0.4 0.6 0.0 0.0 0.0
M3 Ind3 1.0 0.0 0.0 0.0 0.0
M4 Ind3 0.0 0.1 0.9 0.0 0.0
M5 Ind3 0.1 0.9 0.0 0.0 0.0

ft="https://raw.githubusercontent.com/mmollina/SCRI/main/data/B2721_mappoly_prob"
tempfl <- tempfile()
download.file(ft, destfile = tempfl)
dat.prob.mpl <- read_geno_prob(file.in  = tempfl, prob.thres = 0.95, elim.redundant = TRUE)
#> Reading the following data:
#>     Ploidy level: 4
#>     No. individuals:  156
#>     No. markers:  6502
#>     No. informative markers:  6502 (100%)
#>     This dataset contains sequence information.
#>     ...
#>     Done with reading.
#>     Filtering non-conforming markers.
#>     ...
#>     Performing chi-square test.
#>     ...
#>     Done.
unlink(tempfl)

dat.prob.mpl

download.file("https://github.com/mmollina/MAPpoly_vignettes/raw/master/data/BT/sweetpotato_chr1.vcf.gz", destfile = 'chr1.vcf.gz')
dat.dose.vcf = read_vcf(file = 'chr1.vcf.gz', parent.1 = "PARENT1", parent.2 = "PARENT2", verbose = FALSE)
#> Registered S3 method overwritten by 'vegan':
#>   method     from      
#>   rev.hclust dendextend

dat.dose.vcf

library(ggplot2)
dosage_combs = cbind(dat.dose.vcf$dosage.p, dat.dose.vcf$dosage.q)
dc_simplex = apply(dosage_combs,1,function(x) if(all(c(0,1) %in% x) | all(c(dat.dose.vcf$m-1, dat.dose.vcf$m) %in% x)) return(TRUE) else return(FALSE))
dc_dsimplex = apply(dosage_combs,1,function(x) if(all(x == c(1,1)) | all(x == c(dat.dose.vcf$m-1, dat.dose.vcf$m-1))) return(TRUE) else return(FALSE))

dc_simplex[which(dc_simplex == TRUE)] = "simplex"
dc_simplex[which(dc_dsimplex == TRUE)] = 'double simplex'
dc_simplex[which(dc_simplex == 'FALSE')] = 'multiplex'

data_depths = data.frame('Marker depths' = dat.dose.vcf$all.mrk.depth,
                         'Depth classes' = findInterval(dat.dose.vcf$all.mrk.depth, c(200,300,400,500,600,50000)),
                         'Dosage combinations' = dc_simplex, check.names = F)

ggplot(data_depths, aes(fill=`Dosage combinations`, x=`Depth classes`, y=`Marker depths`)) +
  geom_bar(position = 'stack', stat = 'identity') +
  scale_x_continuous(breaks=0:5, labels=c("[0,200)","[200,300)","[300,400)","[400,500)","[500,600)", "> 600"))

## Downloading and reading B2721 fitPoly's probabilistic scores
address <- "https://github.com/mmollina/SCRI/raw/main/data/fitpoly_tetra_call/B2721_scores.zip"
tempfl <- tempfile(pattern = "B2721_", fileext = ".zip")
download.file(url = address, destfile = tempfl)
unzip(tempfl, files = "B2721_scores.dat")
unlink(tempfl)
dat <- read_fitpoly(file.in = "B2721_scores.dat", 
                    ploidy = 4, 
                    parent1 = "Atlantic", 
                    parent2 = "B1829", 
                    verbose = TRUE)
#> Reading the following data:
#>     Ploidy level: 4
#>     No. individuals:  156
#>     No. markers:  7867
#>     No. informative markers:  7093 (90.2%)
#>     ...
#>     Done with reading.
#>     Filtering non-conforming markers.
#>     ...
#>     Performing chi-square test.
#>     ...
#>     Done.

source("https://raw.githubusercontent.com/mmollina/SCRI/main/MAPpoly/get_solcap_snp_pos.R")

# load example dataset from polyRAD
library(polyRAD)
data(exampleRAD_mapping)
exampleRAD_mapping = SetDonorParent(exampleRAD_mapping, "parent1")
exampleRAD_mapping = SetRecurrentParent(exampleRAD_mapping, "parent2")
exampleRAD_mapping = PipelineMapping2Parents(exampleRAD_mapping)
#> Making initial parameter estimates...
#> Generating sampling permutations for allele depth.
#> Updating priors using linkage...
#> Done.

# export to MAPpoly
outfile2 = tempfile()
Export_MAPpoly(exampleRAD_mapping, file = outfile2)

# Read in MAPpoly
mydata_polyrad = read_geno_prob(outfile2)
#> Reading the following data:
#>     Ploidy level: 2
#>     No. individuals:  100
#>     No. markers:  4
#>     No. informative markers:  4 (100%)
#>     This dataset contains sequence information.
#>     ...
#>     Done with reading.
#>     Filtering non-conforming markers.
#>     ...
#>     Performing chi-square test.
#>     ...
#>     Done.

mydata_polyrad

# Load example dataset from updog
library(updog)
data(uitdewilligen)
mout = multidog(refmat = t(uitdewilligen$refmat), 
                sizemat = t(uitdewilligen$sizemat), 
                ploidy = uitdewilligen$ploidy, 
                model = "f1",
                p1_id = colnames(t(uitdewilligen$sizemat))[1],
                p2_id = colnames(t(uitdewilligen$sizemat))[2],
                nc = 4)
mydata_updog = import_from_updog(mout)

mydata_updog

# Downloading VCF files regarding chromosome 1 and 2
download.file("https://github.com/mmollina/MAPpoly_vignettes/raw/master/data/BT/sweetpotato_chr1.vcf.gz", destfile = 'chr1.vcf.gz')
download.file("https://github.com/mmollina/MAPpoly_vignettes/raw/master/data/BT/sweetpotato_chr2.vcf.gz", destfile = 'chr2.vcf.gz')
data1 = read_vcf(file = 'chr1.vcf.gz', parent.1 = "PARENT1", parent.2 = "PARENT2", verbose = FALSE)
#> Registered S3 method overwritten by 'vegan':
#>   method     from      
#>   rev.hclust dendextend
data2 = read_vcf(file = 'chr2.vcf.gz', parent.1 = "PARENT1", parent.2 = "PARENT2", verbose = FALSE)

# See datasets
print(data1)
print(data2)

# Merge datasets
merged_data = merge_datasets(data1, data2)

print(merged_data)

print(dat, detailed = TRUE)

plot(dat)

plot_mrk_info(input.data = dat, mrk = 1979)

plot_mrk_info(input.data = dat, mrk = 'solcap_snp_c2_17752')

print_mrk(dat, mrks = 'solcap_snp_c2_17752')
#> 
#>  solcap_snp_c2_17752
#> ----------------------------------
#>  dosage P1:  2
#>  dosage P2:  2
#> ----
#>  dosage distribution
#>   0   1   2   3   4 mis 
#>   1  29  83  38   5   0 
#> ----
#>  expected polysomic segregation
#>          0          1          2          3          4 
#> 0.02777778 0.22222222 0.50000000 0.22222222 0.02777778 
#> ----------------------------------

3 Filtering and Quality Control

# Filtering dataset by marker
dat <- filter_missing(input.data = dat, type = "marker", 
                               filter.thres = 0.05, inter = FALSE)

print(dat)

# Filtering dataset by individual
dat <- filter_missing(input.data = dat, type = "individual", 
                               filter.thres = 0.025, inter = FALSE)                            

print(dat)

source("https://raw.githubusercontent.com/mmollina/SCRI/main/MAPpoly/get_solcap_snp_pos.R")

pval.bonf <- 0.05/dat$n.mrk
mrks.chi.filt <- filter_segregation(dat, chisq.pval.thres =  pval.bonf, inter = FALSE)
seq.init <- make_seq_mappoly(mrks.chi.filt)

plot(seq.init)

4 Two-point analysis

Once the markers are filtered and selected, we need to compute the pairwise recombination fraction between all of them (two-point analysis). The function est_pairwise_rf is used to estimate all the pairwise recombination fractions between markers in the sequence provided. Since the output object is too big to be fully displayed on the screen, MAPpoly shows a summary. Parallel computation is available and, if you want to use it, you need to define the available number of cores in your machine and also guarantee that you have sufficient RAM memory for that. Remember that it is always good to leave one core available for the system, and be aware that this step will take a while to compute if you have few cores available.

# Defining number of cores
n.cores = parallel::detectCores() - 1

#(~ 9.2 minutes using 23 cores)
all.rf.pairwise <- est_pairwise_rf(input.seq = seq.init, ncpus = n.cores)
#> INFO: Using  23  CPUs for calculation.
#> INFO: Done with 12427605  pairs of markers 
#> INFO: Calculation took: 386.07 seconds

all.rf.pairwise
#>   This is an object of class 'poly.est.two.pts.pairwise'
#>   -----------------------------------------------------
#>   No. markers:                             4986 
#>   No. estimated recombination fractions:   11694917 (94.1%)
#>   -----------------------------------------------------

To assess the recombination fraction between a particular pair of markers, say markers 2204 and 4508, we use the following syntax:

all.rf.pairwise$pairwise$`2204-4508`
#>        LOD_ph         rf       LOD_rf
#> 1-1  0.000000 0.05688641 1.340883e+00
#> 2-1 -2.753497 0.43080101 1.098301e-01
#> 1-2 -2.753497 0.43080101 1.098301e-01
#> 1-0 -2.863464 0.49995672 1.369849e-04
#> 0-1 -2.863464 0.49995672 1.369849e-04
#> 2-2 -4.727155 0.47577396 6.847633e-02
#> 2-0 -4.795632 0.49993129 8.815329e-08
#> 0-2 -4.795632 0.49993129 8.815329e-08
#> 0-0 -4.795880 0.49995600 2.486217e-04
plot(all.rf.pairwise, first.mrk = 2204, second.mrk = 4508)

In this case, 2204-4508 represents the position of the markers in the filtered data set. The name of the rows have the form x-y, where x and y indicate how many homologous chromosomes share the same allelic variant in parents \(P1\) and \(P2\), respectively (see this figure in Mollinari and Garcia (2019) for notation). The first column indicates the LOD Score in relation to the most likely linkage phase configuration. The second column shows the estimated recombination fraction for each configuration, and the third indicates the LOD Score comparing the likelihood under no linkage (\(r = 0.5\)) with the estimated recombination fraction (evidence of linkage).

mat <- rf_list_to_matrix(input.twopt = all.rf.pairwise, ncpus = n.cores)
#> INFO: Using  23  CPUs.
#> INFO: Done with 12427605  pairs of markers 
#> INFO: Operation took: 96.395 seconds

id<-get_genomic_order(seq.init)
s.o <- make_seq_mappoly(id)
plot(mat, ord = s.o$seq.mrk.names, fact = 5)

5 Assembling linkage groups

The function group_mappoly assign markers to linkage groups using the recombination fraction matrix obtained above. The user can provide an expected number of groups or run the interactive version of the function using inter = TRUE. Since in this data set we expect 12 linkage groups (basic chromosome number in potato), we use expected.groups = 12. If the data set provides the chromosome where the markers are located, the function allows comparing the groups obtained using the pairwise recombination fraction and the chromosome information provided using the comp.mat = TRUE.

grs <- group_mappoly(input.mat = mat,
                     expected.groups = 12,
                     comp.mat = TRUE, 
                     inter = FALSE)

grs

Here, we have the 3639 markers distributed in 12 linkage groups. The rows indicate linkage groups obtained using linkage information and the columns are the chromosomes in the reference genome. Notice the diagonal indicating the concordance between the two sources of information. Now, we can plot the resulting marker cluster analysis.

plot(grs)

Once the linkage groups are properly assembled, we use the function make_seq_mappoly to make marker sequences from the group analysis. We will assemble a list with 12 positions, each one containing the corresponding linkage group sequence. Also, we will use only markers allocated in the diagonal of the previous comparison matrix. Thus only markers that were assigned to a particular linkage group using both sources of information will be considered. We will also assemble a smaller two-point object using the functions make_pairs_mappoly and rf_snp_filter to facilitate further parallelization procedures.

# genome correspondence
z <- as.numeric(colnames(grs$seq.vs.grouped.snp)[1:12])
LGS<-vector("list", 12)
for(j in 1:12){
    temp1 <- make_seq_mappoly(grs, j, genomic.info = 1)
    tpt <- make_pairs_mappoly(all.rf.pairwise, input.seq = temp1)
    temp2 <- rf_snp_filter(input.twopt = tpt, diagnostic.plot = FALSE)
    lgtemp <- get_genomic_order(temp2)
    LGS[[z[j]]] <- list(lg = make_seq_mappoly(lgtemp), tpt = tpt)
}

Now, let us print the recombination fraction matrices for each linkage group.

#FIXME:—————– explain that we will construct the map foe one LG and after we will do it for all LGs

6 Estimating the map for a given order

In this section, we will estimate the genetic map for each linkage group. First, we will use the marker order provided by the Solanum tuberosum genome version 4.03. However, when a reference genome is not available or fully colinear to the population’s parents, we need to use de-novo algorithms to order the markers, which will be addressed later in this tutorial.

The estimation of the genetic map for a given order involves the computation of recombination fraction between adjacent markers and finding the linkage phase configuration of those markers in both parents. The core function to perform these tasks in MAPpoly is est_rf_hmm_sequential. This function uses the pairwise recombination fraction as the first source of information to sequentially position allelic variants in specific homologs. For situations where pairwise analysis has limited power, the algorithm relies on the likelihood obtained through a hidden Markov model (HMM) Mollinari and Garcia (2019). Once all markers are positioned, the final map is reconstructed using the HMM multipoint algorithm.

Several arguments are available to control the inclusion and phasing of the markers in the chain. The argument start.set defines the number of initial markers used in an exhaustive search for the most probable configuration. After that, markres are sequentially added to the end of the map. thres.twopt receives the threshold to whether when the linkage phases compared via two-point analysis should be considered, and the HMM analysis should not be used to infer the linkage phase (A. K. A. \(\eta\) in Mollinari and Garcia (2019)). thres.hmm receives the threshold for keeping competing maps computed using HMM (if the two-point analysis was not enough) in the next round of marker insertion. extend.tail indicates the number of markers that should be considered at the end of the chain to compute the multilocus genetic map when inserting a new marker. tol and tol.final receive the desired accuracy to estimate the sub-maps during the sequential phasing procedure and the desired accuracy in the final map. phase.number.limit receives the limit number of linkage phase configurations to be tested using HMM. info.tail is a logical argument: if TRUE it uses the complete informative tail (last markers in the chain that allow all homologous to be distinguished in the parents) of the chain to calculate the likelihood of the linkage phases.

First, as an example, let us estimate the map for linkage group 10. The values used for all arguments were obtained using a balance of processing speed and accuracy of the algorithm. As an exercise, it is interesting to try different values and check out the results. For now, let us stick with the following values (this step can take some hours to finish, depending on chosen parameters):

lg10.map <- est_rf_hmm_sequential(input.seq = LGS[[10]]$lg,
                                start.set = 10,
                                thres.twopt = 10, 
                                thres.hmm = 10,
                                extend.tail = 30,
                                info.tail = TRUE, 
                                twopt = LGS[[10]]$tpt,
                                sub.map.size.diff.limit = 10, 
                                phase.number.limit = 20,
                                reestimate.single.ph.configuration = TRUE,
                                tol = 10e-3,
                                tol.final = 10e-4)
#> ════════════════════════════════════════════════════════════ Initial sequence ══
#> ══════════════════════════════════════════════════ Done with initial sequence ══
#> ══════════════════════════════════ Reestimating final recombination fractions ══
#> ════════════════════════════════════════════════════════════════════════════════

Now, use the functions print and plot to view the results:

print(lg10.map)
#> This is an object of class 'mappoly.map'
#>     Ploidy level:     4 
#>     No. individuals:  144 
#>     No. markers:  210 
#>     No. linkage phases:   1 
#> 
#>     ---------------------------------------------
#>     Number of linkage phase configurations:  1
#>     ---------------------------------------------
#>     Linkage phase configuration:  1
#>        map length:    129.03
#>        log-likelihood:    -1564.04
#>        LOD:       0
#>     ~~~~~~~~~~~~~~~~~~
plot(lg10.map)

Back rectangles indicate the presence of the allelic variant in each one of the four homologous in both parents. One also can print a detailed version of the map and plot a specific map segment using the following code.

print(lg10.map, detailed = TRUE)

plot(lg10.map, left.lim = 30, right.lim = 40, mrk.names = TRUE) 

lg10.map.prob <- est_full_hmm_with_prior_prob(input.map = lg10.map)
plot(lg10.map.prob)

lg10.map.err <- est_full_hmm_with_global_error(input.map = lg10.map, error = 0.05)
plot(lg10.map.err)

lg10.map.updated = update_map(lg10.map.err)
#> Updating map 1
lg10.map.err
#> This is an object of class 'mappoly.map'
#>     Ploidy level:     4 
#>     No. individuals:  144 
#>     No. markers:  210 
#>     No. linkage phases:   1 
#> 
#>     ---------------------------------------------
#>     Number of linkage phase configurations:  1
#>     ---------------------------------------------
#>     Linkage phase configuration:  1
#>        map length:    94.4
#>        log-likelihood:    -2146.67
#>        LOD:       0
#>     ~~~~~~~~~~~~~~~~~~
lg10.map.updated
#> This is an object of class 'mappoly.map'
#>     Ploidy level:     4 
#>     No. individuals:  144 
#>     No. markers:  230 
#>     No. linkage phases:   1 
#> 
#>     ---------------------------------------------
#>     Number of linkage phase configurations:  1
#>     ---------------------------------------------
#>     Linkage phase configuration:  1
#>        map length:    94.4
#>        log-likelihood:    -2146.67
#>        LOD:       0
#>     ~~~~~~~~~~~~~~~~~~

7 Ordering markers using MDS and reestimating the map

So far the map was reestimated using the genomic order. In real situations, unless a genomic information is provided, the markers need to be ordered using an optimization technique. Here, we use the MDS (multidimensional scaling) algorithm, proposed in the context of genetic mapping by Preedy and Hackett (2016). The MDS algorithm requires a recombination fraction matrix, which will be transformed in distances using a mapping function (in this case the Haldane’s mapping function). First, let us gather the pairwise recombination fractions for all three linkage groups:

mt <- lapply(LGS, function(x) rf_list_to_matrix(x$tpt))

Now, for each matrix contained in the object mt, we use the MDS algorithm:

mds.ord <- lapply(mt, mds_mappoly)

Usually at this point, the user can make use of diagnostic plots to remove markers that are disturbing the ordering procedure. Here we didn’t use that procedure, but we encourage the user to check the example in ?mds_mappoly. Now, let us compare the estimated and the genomic orders (feel free to run the last commented line and see interactive plots):

LGS.mds<-vector("list", 12)
for(j in 1:12){
  lgtemp <- make_seq_mappoly(mds.ord[[j]])
  LGS.mds[[j]] <- list(lg = lgtemp, 
      tpt = make_pairs_mappoly(all.rf.pairwise, input.seq = lgtemp))
}

geno.vs.mds <- NULL
for(i in 1:length(LGS.mds)){
  geno.vs.mds<-rbind(geno.vs.mds,
                     data.frame(mrk.names = LGS.mds[[i]]$lg$seq.mrk.names,
                                mds.pos = seq_along(LGS.mds[[i]]$lg$seq.mrk.names),
                                genomic.pos = order(LGS.mds[[i]]$lg$sequence.pos),
                                LG = paste0("LG_", i)))
}

require(ggplot2)
p<-ggplot(geno.vs.mds, aes(genomic.pos, mds.pos)) +
  geom_point(alpha = 1/5, aes(colour = LG)) +
  facet_wrap(~LG) +  xlab("Genome Order") + ylab("Map Order")
p

#plotly::ggplotly(p)

Although several local inconsistencies occurred, the global diagonal pattern indicates a consistent order for all linkage groups using both approaches. Now, let us build the genetic map of linkage group 3 using the MDS order (remember that this can take a while to finish):

lg10.map.mds <- est_rf_hmm_sequential(input.seq = LGS.mds[[10]]$lg,
                                      start.set = 10,
                                      thres.twopt = 10, 
                                      thres.hmm = 10,
                                      extend.tail = 30,
                                      info.tail = TRUE, 
                                      twopt = LGS.mds[[10]]$tpt,
                                      sub.map.size.diff.limit = 10, 
                                      phase.number.limit = 20,
                                      reestimate.single.ph.configuration = TRUE,
                                      tol = 10e-3,
                                      tol.final = 10e-4)
#> ════════════════════════════════════════════════════════════ Initial sequence ══
#> ══════════════════════════════════════════════════ Done with initial sequence ══
#> ══════════════════════════════════ Reestimating final recombination fractions ══
#> ════════════════════════════════════════════════════════════════════════════════

And plot the map:

plot(lg10.map.mds)

It is also possible to compare the maps using both genomic-based and MDS-based orders with the function plot_map_list:

plot_map_list(list(genome = lg10.map, mds = lg10.map.mds), col = c("#E69F00", "#56B4E9"), title = "")

The genomic-based map included the same number of markers but is smaller than the MDS-based map, which indicates a better result. To formally compare the maps, one needs to select the markers that are present in both maps. Interestingly enough, both procedures included the same markers in the final map; however, we provide the code to perform the comparison even if the maps share only a subset of markers:

mrks.in.gen<-intersect(lg10.map$maps[[1]]$seq.num, lg10.map.mds$maps[[1]]$seq.num)
mrks.in.mds<-intersect(lg10.map.mds$maps[[1]]$seq.num, lg10.map$maps[[1]]$seq.num)
if(cor(mrks.in.gen, mrks.in.mds) < 0){
  mrks.in.mds <- rev(mrks.in.mds)
  lg10.map.mds <- rev_map(lg10.map.mds)
}
map.comp.3.gen<-get_submap(input.map = lg10.map, match(mrks.in.gen, lg10.map$maps[[1]]$seq.num), verbose = FALSE)
map.comp.3.mds<-get_submap(input.map = lg10.map.mds, match(mrks.in.mds, lg10.map.mds$maps[[1]]$seq.num), verbose = FALSE)
prob.3.gen<-extract_map(lg10.map)
prob.3.mds<-extract_map(lg10.map.mds)
names(prob.3.gen)<-map.comp.3.gen$maps[[1]]$seq.num
names(prob.3.mds)<-map.comp.3.mds$maps[[1]]$seq.num

matplot(t(data.frame(prob.3.gen,prob.3.mds[names(prob.3.gen)])), 
        type="b", pch="_", col=1, lty=1, lwd = .5, xlab= "", 
        ylab="Marker position (cM)", axes = F)
axis(2)
mtext(text = round(map.comp.3.gen$maps[[1]]$loglike,1), side = 1, adj = 0)
mtext(text = round(map.comp.3.mds$maps[[1]]$loglike,1), side = 1, adj = 1)
mtext(text = "Genomic", side = 3, adj = 0)
mtext(text = "MDS", side = 3, adj = 1)

Please notice that these maps have the same local inversions shown in the dot plots presented earlier. In this case, the log-likelihood of the genomic order is higher than the one obtained using the MDS order. For this linkage group, the genomic-based map was chosen as the best one.

8 Parallel map construction

## ~13.3 min
## Performing parallel computation
my.phase.func<-function(X){
  x<-est_rf_hmm_sequential(input.seq = X$lg,
                                start.set = 10,
                                thres.twopt = 10, 
                                thres.hmm = 10,
                                extend.tail = 30,
                                info.tail = TRUE, 
                                twopt = X$tpt,
                                sub.map.size.diff.limit = 8, 
                                phase.number.limit = 10,
                                reestimate.single.ph.configuration = TRUE,
                                tol = 10e-3,
                                tol.final = 10e-4)
  return(x)
}
cl <- parallel::makeCluster(12)
parallel::clusterEvalQ(cl, require(mappoly))
parallel::clusterExport(cl, "dat")
MAPs <- parallel::parLapply(cl,LGS,my.phase.func)
parallel::stopCluster(cl)

plot_map_list(MAPs, col = "ggstyle")

my.error.func<-function(X){
  x<-est_full_hmm_with_global_error(input.map = X, 
                                    error = 0.05, 
                                    tol = 10e-4, 
                                    verbose = FALSE)
  return(x)
}
cl <- parallel::makeCluster(12)
parallel::clusterEvalQ(cl, require(mappoly))
parallel::clusterExport(cl, "dat")
MAP.err <- parallel::parLapply(cl,MAPs,my.error.func)
parallel::stopCluster(cl)

all.MAPs <- NULL
for(i in 1:12) 
  all.MAPs<-c(all.MAPs, MAPs[i], MAP.err[i])
plot_map_list(map.list = all.MAPs, col = rep(c("#E69F00", "#56B4E9"), 12))

plot_map_list(MAP.err, col = "ggstyle")

plot_genome_vs_map(MAP.err, same.ch.lg = TRUE)

knitr::kable(summary_maps(MAPs))
#> 
#> Your dataset contains removed (redundant) markers. Once finished the maps, remember to add them back with the function 'update_map'.

9 Genotype conditional probabilities

In order to use the genetic map in QTLpoly, one needs to obtain the conditional probability of all possible 36 genotypes along the 12 linkage groups for all individuals in the full-sib population. Let us use the function calc_genoprob_error, which similarly to est_full_hmm_with_global_error allows the inclusion of a global genotyping error:

genoprob.err <- vector("list", 12)
for(i in 1:12)
   genoprob.err[[i]] <- calc_genoprob_error(input.map = MAP.err[[i]], error = 0.05)

Here, a global genotyping error of 5% was used. Each position of any object in the list genoprob.err contains two elements: an array of dimensions \(36 \times number \; of \; markers \times number \; of \; individuals\) and the position of the markers in the maps in centimorgans. Let us display the results for all linkage groups in individual 1:

ind <- 1
op <- par(mfrow = c(3, 4), pty = "s", mar=c(1,1,1,1)) 
for(i in 1:12)
{
  d <- genoprob.err[[i]]$map
  image(t(genoprob.err[[i]]$probs[,,ind]),
        col=RColorBrewer::brewer.pal(n=9 , name = "YlOrRd"),
        axes=FALSE,
        xlab = "Markers",
        ylab = "",
        main = paste("LG", i))
  axis(side = 1, at = d/max(d),
       labels =rep("", length(d)), las=2)
}

par(op)

In this figure, the x-axis represents the genetic map and the y-axis represents the 36 possible genotypes in the full-sib population. The color scale varies from dark purple (high probabilities) to light yellow (low probabilities). With the conditional probabilities computed, it is possible to use the object genoprob.err alongside with phenotypic data as the input of the software QTLpoly, which is an under development software to map multiple QTLs in full-sib families of outcrossing autopolyploid species.

10 Obtaining individual haplotypes

Once ready, the genotypic conditional probabilities can be used to recover any individual haplotype given the map (details described in Mollinari et al. (2020)). To generate this information, one may use the function calc_homoprob to account for the probabilities of each homologous, in all map positions for each individual. For example, let us view the homologous probabilities for chromosome 1 and individual 10:

homoprobs = calc_homoprob(genoprob.err)
#> 
#> Linkage group  1 ...
#> Linkage group  2 ...
#> Linkage group  3 ...
#> Linkage group  4 ...
#> Linkage group  5 ...
#> Linkage group  6 ...
#> Linkage group  7 ...
#> Linkage group  8 ...
#> Linkage group  9 ...
#> Linkage group  10 ...
#> Linkage group  11 ...
#> Linkage group  12 ...

plot(homoprobs, lg = 1, ind = 10)

Using this graphic, it is possible to identify regions of crossing-over occurrence, represented by the inversion of probability magnitudes between homologous from the same parent. It is also possible to view all chromosomes at the same time for any individual by setting the parameter lg = "all". One may use this information to evaluate the quality of the map and repeat some processes with modifications, if necessary.

11 Evaluating the meiotic process

MAPpoly also handles a function to evaluate the meiotic process that guided gamete formation on the studied population. Given genotype conditional probabilities, one may want to account for homologous pairing probabilities and detect the occurrence of preferential pairing, which is possible through the function calc_prefpair_profiles:

prefpairs = calc_prefpair_profiles(genoprob.err)
#> 
#> Linkage group  1 ...
#> Linkage group  2 ...
#> Linkage group  3 ...
#> Linkage group  4 ...
#> Linkage group  5 ...
#> Linkage group  6 ...
#> Linkage group  7 ...
#> Linkage group  8 ...
#> Linkage group  9 ...
#> Linkage group  10 ...
#> Linkage group  11 ...
#> Linkage group  12 ...

The function returns an object of class mappoly.prefpair.profiles, which was saved as prefpairs. This object handles all information necessary to study pairing, such as the probability for each pairing configuration (\(\psi\); see Mollinari and Garcia (2019)) inside each parent. For a more user-friendly visualization of the results, one may want to look at the plot output:

plot(prefpairs)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#save.image("all.analysis.rda")

This graphic shows information about all pairing configurations and their probabilities, the proportion of bivalent/multivalent pairing, and also the p-value for preferential pairing test for all markers inside each parent.